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The low-energy-level macroscopic wave functions of the Bose-Einstein 

condensate(BEC) trapped in a symmetric double-well and a periodic poten- 
tial are obtained by solving the Gross-Pitaevskii equation numerically. The 
ground state tunnel splitting is evaluated in terms of the even and odd wave 
functions corresponding to the global ground and excited states respectively. 
We show that the numerical result is in good agreement with the analytic 
level splitting obtained by means of the periodic instanton method. 

I. INTRODUCTION 

The experimental realization of Bose-Einstein condensation (BEC) in double-well trap 
[1-5] and optical lattice [6-10] has stimulated active research into various aspects of quantum 
tunneling phenomena such as the Josephson junction [11—14], atomic interferometry [4], the 
two- wire waveguide [2], etc. Dalfovo [11] suggested a Josephson like effect by considering a 
confining potential with two wells separated by a barrier. A difference between the chemical 
potentials of the atoms in the two traps can be achieved, for example, by loading a different 
number of atoms in the traps. The first experimental evidence [7] of the oscillating atom 
current was observed instead in an one- dimensional Josephson Junction array realized with 
condensates in a laser standing wave, i.e., an optical lattice [15]. The latest techniques of 
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coherently splitting the condensate by deforming the single optical trap into two wells serve 
as a model system for tunneling in the condensate and provide a perfect demonstration of 
a trapped-atom interferometer [4,5]. 

The coherent tunneling of BEC between double-well traps results in the level splitting of 
the macroscopic ground state and hence the macroscopic coherence, which has been observed 
in interference experiments [1,4]. Recently, the energy-band structure and level splitting due 
to quantum tunneling in two weakly linked condensates in the "phase" representation have 
been evaluated in terms of the periodic instanton method [16], which manifests itself as 
a powerful tool for the calculation of the tunneling rate and a good approximation for the 
dilute boson gas [17]. It, nevertheless, is not able to explicitly take into account the nonlinear 
interaction term in the Gross-Pitaevskii equation (GPE) describing the atom-atom collisions 
in BEC [18]. It is naturally expected that as the number density of atoms in the condensates 
increases, the effect of the nonlinear interaction between atoms would become important. In 
the present paper we solve the GPE numerically in order to have a quantitative evaluation 
of the energy level splitting of the ground state for BEC confined in a symmetric double- well 
trap and an optical lattice and explore the dependence of the energy level splitting on the 
chemical potential and the s-wave scattering length between atoms. 

Although the numerical solution of GPE has been developed into a standard procedure, 
this never prevents us from seeking efficient analytical methods. The advantage of a nonper- 
turbative method is that it gives not only an good description of the tunneling phenomena 
but also a comprehensive physical understanding in the context of quantum field theory. 
The periodic instanton configurations, which have been shown to be a useful tool in several 
areas of research such as spin tunneling, bubble nucleation and string theory, enable also 
the investigation of the finite temperature behavior of these systems. In the case of the 
Bose-Einstein system, it turns out that the periodic instanton method is reliable in evaluat- 
ing the tunnel splitting for BEC trapped in both the double-well and optical lattice in the 
regime of experimental values of the chemical potential and scattering length. The intention 
of this paper is to quantitatively compare the deviation of the instanton result from the 
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exact numerical solution and to address the applicability of the instanton method to actual 
experimental situations. 

We restricted our discussion to a quasi-one- dimensional (Q1D) BEC since it is mathe- 
matically simple, in which the BEC is prepared in optical and magnetic traps by putting 
atoms in cylindrical traps long enough that the one-particle energy-level spacing in the radial 
direction exceeds the interatomic-interaction energy, and the atoms can move effectively in 
the axial direction [19]. 

The paper is organized as follows. In Sec. II we give a brief review of the mean-field 
analysis for BEC trapped in external potentials. In Sec. Ill the energy level splitting is 
derived in terms of the instanton method based on the GPE with, however, the nonlinear 
interaction term included implicitly in the chemical potential. In Sec. IV we present the 
numerical procedure for solving the GPE and evaluation of the ground state level splitting. 
The numerical level splitting is compared with the instanton result in Sec. V and a brief 
summary is given in Sec. VI. 



II. GPE FOR ONE-DIMENSIONAL BOSE GAS 



We are interested in the macroscopic quantum tunneling between the condensates sep- 
arated by potential barriers and the main concern here is how the nonlinear interaction 
between the atoms would affect the level splitting. We begin with the energy functional for 
the condensed bosons of mass m confined in the external potential V(x) given by 
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+ V(x) |^(x)| 2 + ^|^(x)| 4 



(1) 



where the order parameter of the condensate ip (x) are normalized to the number of atoms in 
the condensate / dx \ip (x)\ 2 = N and the ID effective interaction constant Uq = 2h 2 a/ma 2 ± 
[20,21] characterizes the nonlinear interaction in the condensate through s-wave scattering 
length a. Here the radial extension of the ground state wavefunction a_|_ = ^Jh/muj± is a 
typical length scale in the transverse trap with a confinement frequency cu±. The first-order 
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variation of the energy functional leads to the Gross-Pitaevskii equation (GPE), Hip (x) = 
fiip (x), with the chemical potential \i = {ip \H\ ip) /N calculated as the expectation value of 
the Hamiltonian 

The transverse confinement frequency cjj_ should be large compared with /i/h so that the 
condensate is prepared in one dimension. On the other hand, for a less strong transverse 
confinement, atoms will oscillate in all directions which makes the model not exactly solvable. 
As can be seen in the last part of the paper, tunneling would be greatly enhanced in 3D. 

Consider two models of the external potential where atoms can tunnel through the bar- 
riers. The double-well trapping potential of the form 

V dw (x) = V (l-x 2 /xl) 2 (3) 

allows us to investigate the interwell coupling which results in the splitting of the energy 
level. The potential barrier of depth V between the two minima ±x is assumed to be large 
enough so that the overlap between the wave functions relative to the two traps occurs only 
in the classically forbidden region where the interaction can be ignored. The optical lattice 
trapping potential 

V ol {x) = V cos 2 (k x), (4) 

on the other hand, formed by the standing wave laser beams with wavevector k [8], simulates 
the sine-Gordon potential which is widely used in quantum field theory as a periodic field 
model [22] . Quantum tunneling between many wells leads to the formation of energy bands 
due to the spatially periodic potential (4). A path integral calculation [22] was done for 
these quantum tunneling models both for vacuum and excited states neglecting, however, 
the nonlinear interaction. In our previous work [18], the nonlinear interaction between the 
atoms was included in the finite chemical potential and we realized that tunneling occurs 
at the level of chemical potential. Here we will solve the GPE numerically and compare the 
numerical results with the analytical ones. 

4 



For convenience we rescale the energies and distances in units of Tiuq and oscillator 



lengths a = yh/mujo, with uj = yV"(xb)/m being the frequency of small oscillations at 
the bottom of each well Xb in double-well or optical lattice traps. The wave function is 



correspondingly rescaled in units of \JlJa~o so that it remains normalized to N. The GPE 
thus takes the following dimensionless form 



■~ + V(x) + U \f/,(x)\ 2 



if) (x) = nip (x) (5) 



with the potential barrier V and chemical potential /i measured in units of hu and the 
nonlinear interaction parameter becomes U = 2aa Q /a 2 ± . Futhermore, we fix the potential 



parameters as x = v8Vo and k = 1/2 m order to leave us with only one adjustable 
parameter, i.e., the potential depth V . The tunnel splitting depends on these parameters 
xo and ko, which are effectively the separations between the well bottoms, in a similar way 
as on Vq. In the case of an optical lattice, the depth of the barrier is usually measured in 
units of the recoil energy E r = % 2 kl/2m of the atoms. It is, however, not difficult to transfer 
the energy units between this convention and ours. 



III. THE TUNNEL SPLITTING EVALUATED WITH INSTANTON METHOD 

Quantum tunneling between noninteracting particles localized in two adjacent wells with 
macroscopic wave functions ip- leads to an effective energy level splitting A/i, which 
removes the asymptotic degeneracy. The wave functions of the ground state tp e with even 
parity and the first excited state ip Q with odd parity are superpositions of the localized wave 
functions ip + and ■?/>- 

^ e = ty+ + ^-)/V2 (6) 

^ = (^+-^-)/v / 2 (7) 

with energy eigenvalues /i ± A/x/2, respectively. When the interatomic coupling constant 
Uq vanishes, the problem reduces to the solution of a linear Schrondinger equation with the 
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Hamiltonian H — — | J-^ + V (x). The nonlinear interaction increases the chemical potential 
and even for the system at zero temperature, tunneling occurs at a higher level /i than the 
ground state. The tunnel splitting can be found with the instanton method [18,22] and is 
generally expressed as 

A^ = ^ exp [-S(v)] (8) 

7T 

where the imaginary time action S(/j) is calculated through the barrier region (/i < V(x)) 
once from turning point —b to turning point b: 

S(jm) = f fi(V(x)-n)dx (9) 

and the frequency cu(ii) appearing in the prefactor is the frequency of the classical peri- 
odic oscillations at energy fi > V(x) in the classically accessible region with the boundary 
determined by the turning points b and a 

^) = Ta \ ■ ( 10 ) 

The path integral method [23] has been used in the evaluation of the tunneling rate prefactor 
and the barrier V between two wells is assumed to be high enough to safely use the WKB 
wave functions [11] in the calculation of the transition amplitude. For the potential in the 
form of (3) the level splitting reduces to 

W = ^{l + u) 1/2 {S{k)-ulC{k)), (12) 

where /C (k) and S (k) denote the complete Jacobian elliptic integral of the first and second 
kinds respectively. The corresponding parameters are defined as u = ./^, k 2 = and 



k' 2 = 1 — k 2 . When the nonlinear interaction vanishes, the dimensionless chemical potential 
reduces to \i = 1/2 and the above result turns out to be 

^-.,1/2 / 16V \ 
n exp — . 



A/i^^-SVo^exp^-^j. (13) 



which resembles the case of the vacuum instanton. 

Tunneling between many potential wells splits the level further into as many sublevels 
as the number of wells. The Bloch theorem tells us that the eigenvalues of the periodic 
potential exhibits an energy band structure in the tight-binding approximation 

E{9) = cos (6a) (14) 

with A/i the band width of the quantum state with energy \i. Here the Bloch wave vector 6 
is confined to the first Brillouin zone [—71/ a, n/a] of the optical lattice with a lattice constant 
a = A/2. For the periodic potential of the optical lattice (4), the energy band width reads 

^ = ^kwr ( - W) (15) 

W = 4V2V {S{k)-k' 2 IC{k)) (16) 

with k 2 = 1 - h/Vq. 

IV. NUMERICAL PROCEDURE 

There exist various numerical approaches for studying the energy spectrum and dynamics 
of BEC trapped in the external potentials [24]. In the present paper, we solve the GPE 
numerically and find the ground- and first excited-state wave functions ip e (x) , ip Q (x) where 
the corresponding energy expectation values ji e and ji Q are obtained by direct calculation. 
The level splitting A/i = /i e — fi Q is described as a function of parameters iV and Uq. 

We adopt the Gauss-Seidel method to solve eq. (5) numerically and hence start from 
the diffusion equation 

with a diffusion constant of | and a source term p. As t — > +00, the wave function relaxes 
to an equilibrium solution which means that all time derivatives vanish. As a matter of 
fact the diffusion equation eq. (17) is obtained from the NLSE eq. (5) with the time being 
replaced by a negative imaginary time. 



We use the following Crank-Nicholson scheme to discretize eq. (17) by using the space 
step h and time step A 



where ■?/>" = i>(xi,t n ) denotes the exact solution at X; L = ih and t n = nA. The method is 
stable, unitary, and second-order accurate in space and time [25-27]. In a lattice of s points 
eq. (18) represents a tridiagonal set with open boundary conditions or a cyclic tridiagonal 
set with periodic boundary conditions for i — 2, 3, s — 1. For tridiagonal sets, the whole 
solution can be encoded very concisely using the procedures of LU decomposition, forward- 
and back- substitution while for cyclic tridiagonal sets, the procedure of Sherman-Morrison 
Formula is used [25]. For the double- well case we choose the space step h = 0.01, time step 
A = 0.001 and, s = 2400. For the optical lattice, A = 0.02 and s = 2500. The values of 
these parameters are chosen to satisfy the stability criterion of the Crank-Nicholson code. 

We start from the initial, trial wave functions (t = 0) given in eqs.(6) and (7) and choose 
the eigenstates in the non-interacting limit as our trial wave functions such that ifj + , if}~ are 
the degenerate eigenstates in the left- and right-well with the same energy eigenvalue. In our 
procedure, all of the wavefunctions with even parity finally evolve into the lowest eigenstate 
with even parity, i.e., the lower level state for the double well or the bottom of energy band 
for the optical lattice. Similarly those states with odd parity evolve into the lowest state 
with odd parity, i.e., the upper level state for the double well or the top of the energy band 
for the optical lattice. The boundary and normalization conditions are implemented at each 
time step. To test the validity of our code, the numerical wave functions for a stationary 
GPE in a spherical trap is compared with the corresponding results given in [28] and the 
agreement is perfect. 




(18) 



+ (^ +1 -2^ + C-i)l 

-l[Vi(xi)+u m 2 ] 0r +1 +O 



8 



V. NUMERICAL RESULT WITH THE NONLINEAR INTERACTION 



As an example we consider weak-linked condensates of 87 Rb atoms confined in multi- 
traps with frequency uj = 100Hz as in reference [12] and the corresponding oscillator length 
is a = 2.70 x 1CT 4 cm. The s-wave scattering length is in the range 85a,bohr < o, sc < 140ab o hr, 
where cibohr is the Bohr radius [29]. In our analysis we use a sc = lOOa^,,. The transverse 
confinement frequency u± is taken to be 2tt x 250Hz. The corresponding radial extension 
a± = 6.81 x 1CT 4 cm and the interatomic-interaction constant Uo = 0.06 (in units of hujo), 
which corresponds to a weak [30] nonlinear interaction such that we could examine its effect 
on the level splitting. We always measured energies in units of huo and lengths in units of 
the oscillator length a , so V is all we need. 

As a comparison we first of all deal with the "noninteracting" case. When V = 5, the 
analytical tunnel splitting of the instanton approach given in eq. (11) is A/i = 3.74 x 10~ n 
and our numerical result is A/i = 3.60 x 10~ n . This again proves the validity of our 
numerical simulation. The corresponding wave function is shown in Fig. la. In this paper 
we are mainly interested in the ground state, and for this purpose the choice of s = 2400 
is seen to be adequate for most of the calculations. In Fig. lb we show the profiles of the 
even (solid line)- and odd (dotted line) wave functions ip ej0 (x) for 344 atoms confined in the 
trap with height V = 5, which are the even- and odd- eigenstates of the Hamiltonian with 
the nonlinear interaction term Uq \ip {x)\ 2 . Also the Thomas- Fermi approximation ipTF = 
[(/i — V (x)) /Uq] 1 ^ 2 for V (x) < /i is given by the dashed lines. 

Now we turn to examine the difference of level splitting obtained from the instanton 
method and the numerical simulation for the double-well trap. In Fig. 2, A/i on a logarithmic 
scale to base 10 as a function of the chemical potential is depicted for barrier heights, V = 4 
and Vo = 5 . Results from both the numerical simulation (solid line) and instanton approach 
(dotted line) exhibit an enhancement of the tunneling with increasing chemical potential. To 
show how large the difference is, in the insets, we plot the splitting divided by the exponential 
factor of the analytical result. The results of the two approaches always have the same order 
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of magnitude and are close to each other. The periodic instanton method evaluates fairly 
well the tunneling splitting even if the nonlinear interaction is included. Quantitatively, 
however, it always over-estimates the splitting for the whole range of the chemical potential. 
These can be seen more clearly from the interacting constant U dependence of the level 
splitting displayed in Fig. 3 and its inset for V = 5. We emphasize here that a peculiar 
feature of this periodic instanton result for the level splitting is that the prefactor depends 
on the chemical potential as displayed in the insets of Fig. 2. This result is important, as 
has been shown in [18]. 

For the case of the optical lattice potential we still choose the same parameters as in the 
case of the double- well. The nonlinear interaction constant between atoms in the same well, 
Uq = 0.06 for the repulsive interaction. The typical profile of the condensate wave function 
if) e (x), shown by the solid line, and ip (x), shown by the dotted lines, are plotted in Fig. 
4 for Vo = 5 and 150 atoms in each well. The even wave function is symmetric, and the 
odd wave function is maximally antisymmetric, i.e., the wavefunction segment in each well 
is antisymmetric with respect to those of its neighbors. Fig. 5a and 5b display the chemical 
potential dependence of the level splitting for Vo = 5 and Vo = 4, as determined by the 
method outlined above. Again, we find the results from the periodic instanton method and 
GPE have the same order of magnitude. 

By comparing Fig. 2 and Fig. 5, it is shown that when the nonlinear interaction between 
atoms is included, the level splitting is smaller than the instanton result in the double-well 
case, but the splitting is larger in the optical lattice case. This distinction depends on the 
structure of the trapping potentials, for example, atoms in the optical lattice can tunnel 
through the barrier in two directions, while the tunnel path for those in the double-well is 
one-way only. This makes the periodic instanton result different from the double well case. 
In spite of this, the periodic instanton method remains good enough to describe the level 
splitting for the BEC. 

Recently a single bosonic Josephson junction [31] has been implemented by two weakly 
linked BEC in a double-well potential. In previously reported realizations of condensates in 
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double-well potentials [5] the time scale of tunneling dynamics was in the range of thousands 
of seconds. In contrast, their setup allows the realization of nonlinear tunneling times on 
the order of 50ms, which makes the direct observation of the nonlinear dynamics in a single 
bosonic Josephson junction possible for the first time. We emphasize here the distinction 
between tunneling in 3D and that in quasi-lD systems. The important parameter, the 
"tunneling matrix element" K [12] between two condensates is related to the energy level 
splitting through K = A/x/2. The period of oscillation T can be obtained by numerically 
integrating eq. (1) in Ref. [31] for each given parameter A = NU /2K. In the 3D case, K 
is often assumed to be of the order O.lnK or 25Hz, which gives a relative small value for 
A ~ 10. In contrast, in the ID case, the tight confinement in the other 2 directions would 
suppress drastically the tunneling and make the link between the condensates even weaker. 
According to our calculation the parameter K obviously takes typical values of 10~ 3 ~ 1Hz 
and A may be as large as 10 3 ~ 10 6 . As a consequence, the atoms tend to be trapped in the 
potential wells and the observation of Josephson oscillation becomes almost impossible, e.g., 
the initial population imbalance must be less than 0.06 for A ~ 10 3 . A simple calculation 
shows that for N = 600, U = 0.06, V = 5, the Josephson oscillation may be observed with 
a period 2ms, which is less than that in the 3D experiments. 

According to our calculation, the energy splitting is very small compared with the chem- 
ical potential of the system and increases exponentially with the chemical potential. The 
tunneling effect gives rise to the macroscopic phase coherence of BEC across the barriers, 
which results in the significant observable interference phenomena between different BEC 
segments. 

VI. SUMMARY 

In this paper we have investigated the tunnel splitting of the ground state for a weakly- 
linked BEC trapped in double-well potential and in optical lattice by solving the Gross- 
Pitaevskii equations numerically. It turns out that the periodic instanton method is a 
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reliable tool in the evaluation of the tunnel splitting for the BEC. The reason for this is 
that for the quantum tunneling through the potential barrier, the nonlinear interaction is 
negligibly small and contributes overall to a finite chemical potential. Our numerical scheme 
could easily be improved for the investigation of the dynamical behavior of the condensates 
in multi-well potentials. 
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Figure Captions: 

1. Level splitting and wavefunctions of condensates confined in a double well potential for 
Vo = 5. Panel (a) The symmetric (full line) and antisymmetric (dotted lines) wave functions 
of the ground state with U = (noninteracting case). Panel (b) The same as Panel (a) 
but for Uq = 0.06. The Thomas-Fermi solution ip TF = [(fi — V (x)) /Uq] 1 ^ 2 is shown by the 
dashed lines. 

2. Level splitting as a function of chemical potential /x for two values of the barrier 
height V — 5 (a) and V = 4 (b). Insets: The splitting divided by the exponential factor of 
the analytical result in the insets. Solid lines: GPE results Afi GPE / exp(— W), dotted lines: 
periodic instanton results A// Instanton / exp(— W). 

3. The interaction constant Uq dependence of level splitting for condensates in a double- 
well for Vo = 5. The inset shows the splitting divided by the exponential factor of analytical 
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result. Solid lines: GPE results, A/i G pe/ exp(— W), dotted lines: periodic instanton results 

A/Xl ns t a nton/ 6Xp(-W). 

4. The symmetric (full line) and antisymmetric (dotted lines) wave functions in an 
optical lattice for V = 5. 

5. Level splitting as a function of chemical potential /x in an optical lattice for V — 5 
(a) and Vq = 4 (b). The insets show the splitting divided by the exponential factor of 
the analytical result. Solid lines: GPE results Afi GPE / exp(— W), dotted lines: periodic 
instanton results A/ii nstant0 n/ exp(— W). 



15 




Fig.1 



8 ,- 
6 - 
4 - 



- 
-2 - 



i 

1 1 



'J 



I". 



I' .1 



U >J l 1 <j 



I 1 -J l' y 



nil 



ill 



U'J 



i' j 



even state 
odd state 



V 'J 



."J 



-6 - 



-300 



-200 



-100 







100 



200 



300 



